/*
*  W. Elliott
*  W. T. Rankin
*
*  Copyright (c) 1994 Duke University
*     All rights reserved.				
*/

#ifndef _multipole_h
#define _multipole_h

#include "mpole_misc.h"
#include "vektor.h"
#include "particle.h"
#include "dpmta_legendre.h"

/* global constants */
#define SMALL_THETA 1.0e-10
#define M_PI_2      1.57079632679489661923

extern int Mtermp;
extern int particleLEAF;	//number of particle in a LEAF 

/* global data structures */
static Mtype   Y_C;
static Mtype   Hm2l;
static Mtype   L;
static Kompleks *Yxy;
static Real    **Fcoeff;
static Real    **Gcoeff;
static Real    **LegPoly;
static int     *PFindex;


/****************************************************************
*
*  this procedure performs the one-time mpe initialization of
*  global constants and arrays
*
*/
int MultipoleSetup( int Mp );

/*****************************************************************
*
*  AddMultipoleC calculates the Coulomb multipole expansion of a particle 
*  of charge q, at location <vect> relative to the origin. The result is
*  accumulated into the blank multipole expansion array M. The array size
*  is p.
*
*/
void AddMultipoleC(Kompleks **M, int p, Real q, vektor vect);

/****************************************************************
*
*  Calc_M2M() - shift the multipole expansion for a child cell and add
*    it to the multipole expansion of the parent cell.
*
*  w.t.rankin
*
*/

int Calc_M2M( Node* ccell, Node* pcell);

/*****************************************************************
*
*  M2M_C() translates a Coulomb multipole expansion M1 to another 
*  expansion M2, both size p, along the vector <vect>
*
*/

void M2M_C( Kompleks **M1, Kompleks **M2 , int p , vektor vect );

/****************************************************************
*
*  Calc_M2L() - shift the multipole expansion for a remote cell and add
*    it to the local expansion of the local cell.
*
*/
int Calc_M2L(Node* rcell , Node* lcell );

/*****************************************************************
*
*  M2L_C converts a Coulomb multipole expansion M in to a local 
*  expansion L, shifting the original expansion along the vector <vect>.
*  Both expansions are of size p.
*
*/

int M2L_C( Kompleks **M, Kompleks **L , int p , vektor vect );

/****************************************************************
*
*  Calc_L2L() - take the local expansion for a parent cell and
*     add it to the local expansion of a child cell.
*/

void Calc_L2L( Node* pcell, Node* ccell);

void makeG( int , vektor );
void makeF( int , vektor );


/*****************************************************************
*
*  L2L_C translates a local expansion L1 to another local expansion L2, 
*  both of size p, along vector <vect>
*
*/

int L2L_C( Kompleks **L1, Kompleks**L2 , int p , vektor vect ) ;


/*****************************************************************
*
*  Force_C computes the Coulomb force and potential on the particle with
*  charge q at position <vect>, due to the local expansion L of size p. 
*  The results are returned as a potential 'rnt_pot', and as a vector
*  force 'rnt_force'.
*
*/
void Force_C( Kompleks **L, int p, double q, vektor vect, double *rtn_pot, vektor *rtn_force);

/*****************************************************************
*
*  allocate and initialize static data structures needed by the
*  multipole calculations.
*
*/
void Cinit(int p );

/*****************************************************************
*
*  allocates the multipole expansion of size p leaving the resulting
*  matrix in Mptr.
*
*/

int Calloc(Kompleks ***Mptr , int p ) ;

/****************************************************************
*
*  CMclear(M,p) -
*     clears the multipole expansion of size p pointer to by Mptr
*
*/

int CMclear( Kompleks **M , int p );

/****************************************************************
*
*  CMsum(M1,M2,p) -
*     adds the multipole expansion M1 (of size p) to M2
*
*/

int CMsum(Kompleks **M1,Kompleks **M2,int p);

static void Fourier_C( int p, Real b);


void Clear_mpole( Node* thiscell);

void Clear_local( Node* thiscell);


void display_expansion ( Kompleks **M , int p );

void Display_mpole( Node* thiscell );

void Display_local( Node* thiscell );

#endif
